"""
Main module for differential operators across all coordinate systems.

This module provides a unified interface for creating and using differential operators
in PDEJulia.jl. It automatically dispatches to the appropriate coordinate system
implementation and handles optimization, caching, and boundary conditions.

Usage:
    # Create operators
    laplacian = create_operator(:laplacian, grid, boundary_conditions)
    gradient = create_operator(:gradient, grid, boundary_conditions, dim=1)
    
    # Apply operators
    apply_operator!(laplacian, input_field, output_field)
    
    # High-level interface
    result = laplace(field, grid, bcs)

Supported operators:
- :laplacian - Scalar Laplacian ∇²u
- :gradient - Gradient ∇u (specify dimension)
- :divergence - Divergence ∇·v for vector fields
- :curl - Curl ∇×v for vector fields  
- :vector_laplacian - Vector Laplacian ∇²v
- :biharmonic - Biharmonic ∇⁴u
- :poisson_solver - Poisson equation solver ∇²u = f

Supported coordinate systems:
- CartesianGrid: Standard Cartesian coordinates
- SphericalSymGrid: Spherically symmetric coordinates
- CylindricalSymGrid: Cylindrically symmetric coordinates  
- PolarSymGrid: Polar symmetric coordinates
"""

# Include all operator modules
include("common.jl")
include("cartesian.jl")
include("spherical_sym.jl")
include("cylindrical_sym.jl")
include("polar_sym.jl")
include("boundary_conditions.jl")
include("optimization.jl")

using LinearAlgebra

# Export main interface functions
export create_operator, apply_operator!, laplace, gradient, divergence, curl
export poisson_solve, create_operator_set, OperatorSet
export get_available_operators, validate_operator_request

"""
    OperatorSet

Container for a complete set of operators for a given grid and boundary conditions.
"""
struct OperatorSet
    grid::AbstractGrid
    boundary_conditions::BoundaryConditionSet
    operators::Dict{Symbol, Function}
    optimization_level::Symbol
    
    function OperatorSet(grid::AbstractGrid, bcs::BoundaryConditionSet; 
                        optimization_level::Symbol = :medium)
        new(grid, bcs, Dict{Symbol, Function}(), optimization_level)
    end
end

"""
    create_operator(operator_type::Symbol, grid::AbstractGrid, bcs; kwargs...)

Create a differential operator for the specified grid and boundary conditions.

Args:
    operator_type: Type of operator (:laplacian, :gradient, :divergence, etc.)
    grid: Grid on which to define the operator
    bcs: Boundary conditions
    kwargs: Additional operator-specific parameters

Returns:
    Function that applies the operator
"""
function create_operator(operator_type::Symbol, grid::AbstractGrid, bcs; kwargs...)
    # Validate inputs
    validate_operator_request(operator_type, grid, bcs; kwargs...)
    
    # Check for cached operator first
    cached_op = get_cached_operator(grid, bcs, operator_type)
    if cached_op !== nothing
        return cached_op.operator
    end
    
    # Create new operator based on grid type and operator type
    grid_type = typeof(grid)
    
    if grid_type <: CartesianGrid
        operator = create_cartesian_operator(operator_type, grid, bcs; kwargs...)
    elseif grid_type <: SphericalSymGrid
        operator = create_spherical_operator(operator_type, grid, bcs; kwargs...)
    elseif grid_type <: CylindricalSymGrid
        operator = create_cylindrical_operator(operator_type, grid, bcs; kwargs...)
    elseif grid_type <: PolarSymGrid
        operator = create_polar_operator(operator_type, grid, bcs; kwargs...)
    else
        error("Unsupported grid type: $grid_type")
    end
    
    # Apply optimizations
    optimization_level = get(kwargs, :optimization_level, :medium)
    optimized_operator = optimize_operator(operator, grid, optimization_level)
    
    # Cache the operator
    cache_operator(() -> optimized_operator, grid, bcs, operator_type)
    
    return optimized_operator
end

"""
    create_cartesian_operator(operator_type, grid, bcs; kwargs...)

Create operator for Cartesian grids.
"""
function create_cartesian_operator(operator_type::Symbol, grid::CartesianGrid, bcs; kwargs...)
    if operator_type == :laplacian
        conservative = get(kwargs, :conservative, false)
        order = get(kwargs, :order, 2)
        return make_laplace(grid; conservative=conservative, order=order)
        
    elseif operator_type == :gradient
        dim = get(kwargs, :dim, 1)
        order = get(kwargs, :order, 2)
        return make_gradient(grid, dim; order=order)
        
    elseif operator_type == :divergence
        return make_divergence(grid)
        
    elseif operator_type == :curl
        return make_curl(grid)
        
    elseif operator_type == :vector_laplacian
        return make_vector_laplace(grid)
        
    elseif operator_type == :biharmonic
        return make_biharmonic(grid)
        
    elseif operator_type == :poisson_solver
        return make_poisson_solver(grid, bcs)
        
    elseif operator_type == :fft_laplacian
        return fft_laplace(grid)
        
    else
        error("Unsupported Cartesian operator: $operator_type")
    end
end

"""
    create_spherical_operator(operator_type, grid, bcs; kwargs...)

Create operator for spherical grids.
"""
function create_spherical_operator(operator_type::Symbol, grid::SphericalSymGrid, bcs; kwargs...)
    if operator_type == :laplacian
        conservative = get(kwargs, :conservative, false)
        return make_laplace_spherical(grid; conservative=conservative)
        
    elseif operator_type == :gradient
        method = get(kwargs, :method, "central")
        return make_gradient_spherical(grid; method=method)
        
    elseif operator_type == :divergence
        return make_divergence_spherical(grid)
        
    elseif operator_type == :vector_laplacian
        return make_vector_laplace_spherical(grid)
        
    elseif operator_type == :poisson_solver
        return make_poisson_solver_spherical(grid, bcs)
        
    else
        error("Unsupported spherical operator: $operator_type")
    end
end

"""
    create_cylindrical_operator(operator_type, grid, bcs; kwargs...)

Create operator for cylindrical grids.
"""
function create_cylindrical_operator(operator_type::Symbol, grid::CylindricalSymGrid, bcs; kwargs...)
    if operator_type == :laplacian
        conservative = get(kwargs, :conservative, false)
        return make_laplace_cylindrical(grid; conservative=conservative)
        
    elseif operator_type == :gradient
        dim = get(kwargs, :dim, 1)
        method = get(kwargs, :method, "central")
        return make_gradient_cylindrical(grid, dim; method=method)
        
    elseif operator_type == :divergence
        return make_divergence_cylindrical(grid)
        
    elseif operator_type == :vector_laplacian
        return make_vector_laplace_cylindrical(grid)
        
    elseif operator_type == :poisson_solver
        return make_poisson_solver_cylindrical(grid, bcs)
        
    else
        error("Unsupported cylindrical operator: $operator_type")
    end
end

"""
    create_polar_operator(operator_type, grid, bcs; kwargs...)

Create operator for polar grids.
"""
function create_polar_operator(operator_type::Symbol, grid::PolarSymGrid, bcs; kwargs...)
    if operator_type == :laplacian
        conservative = get(kwargs, :conservative, false)
        return make_laplace_polar(grid; conservative=conservative)
        
    elseif operator_type == :gradient
        method = get(kwargs, :method, "central")
        return make_gradient_polar(grid; method=method)
        
    elseif operator_type == :divergence
        return make_divergence_polar(grid)
        
    elseif operator_type == :vector_laplacian
        return make_vector_laplace_polar(grid)
        
    elseif operator_type == :poisson_solver
        return make_poisson_solver_polar(grid, bcs)
        
    else
        error("Unsupported polar operator: $operator_type")
    end
end

"""
    apply_operator!(operator_func::Function, input_field, output_field)

Apply operator to input field and store result in output field.
"""
function apply_operator!(operator_func::Function, input_field, output_field)
    return operator_func(input_field, output_field)
end

"""
    laplace(field, grid, bcs; kwargs...) -> Array

High-level interface for Laplacian operator.
"""
function laplace(field, grid, bcs; kwargs...)
    laplacian = create_operator(:laplacian, grid, bcs; kwargs...)
    result = similar(field)
    laplacian(field, result)
    return result
end

"""
    gradient(field, grid, bcs, dim::Int; kwargs...) -> Array

High-level interface for gradient operator.
"""
function gradient(field, grid, bcs, dim::Int; kwargs...)
    grad_op = create_operator(:gradient, grid, bcs; dim=dim, kwargs...)
    result = similar(field)
    grad_op(field, result)
    return result
end

"""
    divergence(vector_field, grid, bcs; kwargs...) -> Array

High-level interface for divergence operator.
"""
function divergence(vector_field, grid, bcs; kwargs...)
    div_op = create_operator(:divergence, grid, bcs; kwargs...)
    result_shape = size(vector_field)[1:end-1]  # Remove vector component dimension
    result = zeros(result_shape...)
    div_op(vector_field, result)
    return result
end

"""
    curl(vector_field, grid, bcs; kwargs...) -> Array

High-level interface for curl operator.
"""
function curl(vector_field, grid, bcs; kwargs...)
    curl_op = create_operator(:curl, grid, bcs; kwargs...)
    
    # Determine output shape based on input dimensions
    if ndim(grid) == 2
        result_shape = size(vector_field)[1:end-1]  # Scalar curl in 2D
        result = zeros(result_shape...)
    else  # 3D
        result_shape = size(vector_field)  # Vector curl in 3D
        result = zeros(result_shape...)
    end
    
    curl_op(vector_field, result)
    return result
end

"""
    poisson_solve(rhs, grid, bcs; kwargs...) -> Array

High-level interface for Poisson equation solver: ∇²u = rhs
"""
function poisson_solve(rhs, grid, bcs; kwargs...)
    solver = create_operator(:poisson_solver, grid, bcs; kwargs...)
    result = similar(rhs)
    solver(rhs, result)
    return result
end

"""
    create_operator_set(grid, bcs; operators=:all, kwargs...) -> OperatorSet

Create a complete set of operators for a grid.
"""
function create_operator_set(grid, bcs; operators=:all, kwargs...)
    op_set = OperatorSet(grid, bcs; kwargs...)
    
    if operators == :all
        # Create all supported operators for this grid type
        available_ops = get_available_operators(typeof(grid))
    elseif isa(operators, Vector{Symbol})
        available_ops = operators
    else
        error("Invalid operators specification: $operators")
    end
    
    for op_type in available_ops
        try
            op_set.operators[op_type] = create_operator(op_type, grid, bcs; kwargs...)
            println("Created $op_type operator")
        catch e
            println("Warning: Could not create $op_type operator: $e")
        end
    end
    
    return op_set
end

"""
    get_available_operators(grid_type::Type) -> Vector{Symbol}

Get list of operators available for a given grid type.
"""
function get_available_operators(grid_type::Type)
    if grid_type <: CartesianGrid
        return [:laplacian, :gradient, :divergence, :curl, :vector_laplacian, 
                :biharmonic, :poisson_solver, :fft_laplacian]
    elseif grid_type <: SphericalSymGrid
        return [:laplacian, :gradient, :divergence, :vector_laplacian, :poisson_solver]
    elseif grid_type <: CylindricalSymGrid
        return [:laplacian, :gradient, :divergence, :vector_laplacian, :poisson_solver]
    elseif grid_type <: PolarSymGrid
        return [:laplacian, :gradient, :divergence, :vector_laplacian, :poisson_solver]
    else
        return Symbol[]
    end
end

"""
    validate_operator_request(operator_type, grid, bcs; kwargs...)

Validate that the operator request is valid for the given grid and boundary conditions.
"""
function validate_operator_request(operator_type::Symbol, grid, bcs; kwargs...)
    # Check if operator is supported for this grid type
    available_ops = get_available_operators(typeof(grid))
    if !(operator_type in available_ops)
        error("Operator $operator_type not supported for grid type $(typeof(grid))")
    end
    
    # Validate boundary condition compatibility
    validate_bc_compatibility(bcs, operator_type)
    
    # Check operator-specific requirements
    if operator_type == :gradient
        dim = get(kwargs, :dim, 1)
        if dim < 1 || dim > ndim(grid)
            error("Invalid gradient dimension $dim for $(ndim(grid))D grid")
        end
    end
    
    if operator_type == :fft_laplacian
        # Check for periodic boundary conditions
        if !all(bc -> isa(bc, PeriodicBC), bcs.conditions)
            error("FFT Laplacian requires periodic boundary conditions")
        end
    end
    
    return true
end

"""
    Base.show(io::IO, op_set::OperatorSet)

Pretty printing for operator sets.
"""
function Base.show(io::IO, op_set::OperatorSet)
    println(io, "OperatorSet for $(typeof(op_set.grid))")
    println(io, "  Grid size: $(size(op_set.grid))")
    println(io, "  Optimization level: $(op_set.optimization_level)")
    println(io, "  Available operators:")
    for (name, _) in op_set.operators
        println(io, "    • $name")
    end
end

"""
    demonstrate_operators(grid_type::Symbol = :cartesian)

Demonstrate operator usage with example calculations.
"""
function demonstrate_operators(grid_type::Symbol = :cartesian)
    println("Demonstrating operators for $grid_type grid...")
    
    # Create appropriate grid
    if grid_type == :cartesian
        grid = CartesianGrid((-1.0, 1.0), (-1.0, 1.0), (50, 50))
        bcs = create_bc_set_cartesian_2d(
            DirichletBC(0.0), DirichletBC(0.0), 
            DirichletBC(0.0), DirichletBC(0.0)
        )
    elseif grid_type == :spherical
        grid = SphericalSymGrid(50, (0.1, 2.0))
        bcs = create_bc_set_spherical_1d(NeumannBC(0.0), DirichletBC(0.0))
    else
        error("Unsupported grid type for demonstration: $grid_type")
    end
    
    # Create test field
    if grid_type == :cartesian
        x_coords, y_coords = coordinate_arrays(grid)
        test_field = [exp(-(x^2 + y^2)) for x in x_coords, y in y_coords]
        println("  Test field: 2D Gaussian")
        
        # Demonstrate operators
        laplacian_result = laplace(test_field, grid, bcs)
        println("  Laplacian computed: min=$(minimum(laplacian_result)), max=$(maximum(laplacian_result))")
        
        grad_x = gradient(test_field, grid, bcs, 1)
        grad_y = gradient(test_field, grid, bcs, 2)
        println("  Gradient computed: x-component range $(extrema(grad_x))")
        println("                     y-component range $(extrema(grad_y))")
        
    elseif grid_type == :spherical
        r_coords = coordinate_arrays(grid)[1]
        test_field = exp.(-r_coords.^2)
        println("  Test field: 1D radial Gaussian")
        
        laplacian_result = laplace(test_field, grid, bcs)
        println("  Spherical Laplacian: min=$(minimum(laplacian_result)), max=$(maximum(laplacian_result))")
        
        grad_r = gradient(test_field, grid, bcs, 1)
        println("  Radial gradient: range $(extrema(grad_r))")
    end
    
    println("Operator demonstration complete!\n")
end